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1.  Summary 


The  FATEPEN  model  prediets  the  penetration  of  a  mass  striking  a  target  plate  for 
a  variety  of  shapes,  ineluding  eylinders  and  euboids — among  others.  Crueial  to 
the  use  of  the  model  is  a  good  estimate  of  not  only  the  mass  and  veloeity  but  also 
the  impaet  orientation  in  terms  of  piteh,  yaw,  and  roll.  Yaw  eards  and  orthogonal 
X-rays  ean  provide  estimates  of  the  impaet  projeeted  area,  but  the  model  makes  use 
of  the  impaet  angle,  whieh  is  defined  to  be  the  minimum  angle  that  a  penetrator 
faee  makes  with  the  target  plate.  This  report  addresses  this  issue  by  ealeulating  the 
piteh-yaw-roll  rotation  sequenee  that  will  bring  about  the  orientation  at  impaet  from 
orthogonal  projeeted  area  measurements.  It  is  shown  that  the  impaet  angle  is  uniquely 
determined  in  the  ease  of  cylinders,  but  that  is  not  the  case  for  cuboids.  Furthermore, 
for  cuboids  multiple  impact  angles  for  the  same  projected  area  measurements  can 
lead  to  significantly  different  FATEPEN  predictions.  This  leads  to  the  conclusion 
that  cylinders  and  not  cuboids  should  be  used  for  EATEPEN  validation. 

2.  Introduction 

The  EATEPEN^  model  is  a  Eortran  code  for  simulating  penetration  of  fragments, 
long  rods,  and  projectiles  into  target  plates  and  predicting  penetration,  perforation, 
residual  mass,  and  residual  velocity.  It  provides  the  following  shapes:  sphere,  right 
circular  cylinder  (RCC),  round-nose  cylinder,  sharp-nose  cylinder,  tapered/truncated 
cylinder,  and  rectangular  parallelepiped  (RPP).  This  report  is  limited  to  RCCs,  which 
we  will  simply  call  cylinders,  and  RPPs,  which  we  will  call  cuboids.  Unlike  the 
THOR^  model,  which  only  requires  an  impact  presented  area,  the  EATEPEN  model 
requires  a  specific  shape  and  orientation. 

The  impact  presented  area,  or  projected  area,  is  a  key  parameter  that  affects  the 
penetration  process  since  the  penetration  depth  scales  with  the  mass  per  unit  pre¬ 
sented  area  for  a  given  striking  velocity.  But  the  EATEPEN  model  goes  beyond 
the  presented  area  and  also  requires  the  angle  that  each  flat  face  of  the  striking 
mass  makes  with  the  target  plate.  This  means  that  we  need  to  account  for  both  the 
presented  area  and  the  orientation  of  the  fragment  to  use  EATEPEN  properly. 
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3.  Methods,  Assumptions,  and  Procedures 


We  first  define  a  standard  orientation  of  the  eylinder  or  euboid  as  the  starting 
orientation  and  then  deseribe  the  operational  proeedure  to  give  it  a  speeifie  final 
orientation.  This  standard  orientation  will  be  with  the  eenter  at  the  origin  of  a  right- 
handed  Cartesian  {x,  y,  z)  eoordinate  system  and  the  length  of  the  eylinder  or  euboid 
aligned  with  the  z-axis.  (The  target  plate  is  taken  to  be  parallel  to  the  x-y  plane,  and 
the  fragment  veloeity  veetor  is  along  the  negative  2;  axis.)  The  operational  procedure 
will  be  in  terms  pitch,  yaw,  and  roll,  where  pitch  is  a  counterclockwise  rotation 
about  the  x-axis,  yaw  is  a  counterclockwise  rotation  about  the  y-axis,  and  roll  is  a 
counterclockwise  rotation  about  the  2;-axis.*  For  convenience,  we  also  define  i,  j, 
and  k  to  be  unit  vectors  along  the  x,  y,  and  z  axis,  respectively,  as  shown  in  Fig.  1. 

y 


Fig.  1.  Description  of  pitch,  yaw,  and  roll 


Pitch,  yaw,  and  roll  are  then  described  by  the  operators  i?i(0p),  and 

where  the  unit  vector  subscript  denotes  the  axis  of  rotation,  and  (j)p,  (j)y,  and  0^  are 
pitch,  yaw,  and  roll  angles,  respectively.  A  rotation  sequence  accounts  for  the  fact 
that  the  unit  vectors  are  also  transformed  by  prior  rotations.  Thus,  a  pitch-yaw-roll 
rotation  sequence  is  described  by 

R  =  R^„{(l)r)Ry{(t)y)R\{(t)p),  (1) 

where  the  operators  are  applied  from  right  to  left — first  pitch,  then  yaw,  then  roll — 
and  where  j'  is  the  transformed  j  vector  after  the  pitch  rotation,  and  k"  is  the  doubly 
transformed  k  vector  after  the  pitch  and  yaw  rotations.  Now  we  make  use  of  the 
fact  that  we  generate  the  same  net  rotation  if  we  apply  the  sequence  in  reverse  order 

*To  avoid  any  confusion,  note  that  FATEPEN  uses  a  pitch-yaw-roll  rotation  sequence,  which  we 
follow  in  this  report,  but  much  of  the  aerospace  industry  uses  a  yaw-pitch-roll  rotation  sequence. 
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about  the  fixed  axes.*  This  means  that 


R  =  Ry,{(f)r)Ry{(l>y)Ri{(l>p)  =  Ri{(i>p)Rfi(l)y)Ri,{(l>r) 


(2) 


and  we  only  have  to  concern  ourselves  with  rotations  about  fixed  axes.  Table  1  lists 
all  3  X  3  =  9  combinations  of  rotations  applied  to  the  unit  vectors. 


Table  1 .  Pitch,  yaw,  and  roll  rotations  of  unit  vectors 


Rotation 

i 

J 

k 

Pitch;  Ri{(j)p) 

i 

cos  (ppj  +  sin  (ppk 

—  sin  ppj  +  cos  Ppk 

Yaw:  Rfi(j>y) 

cos  (pyl  —  sin0yk 

J 

sin  pyi  +  cos  pyk 

Roll: 

cos  (pR  +  sin  prj 

—  sin  PR  -f  cos  prj 

k 

3.1  Orientation  of  a  Cylinder 

The  cylinder  has  a  length  L  and  diameter  D.  In  standard  orientation  the  center  of  the 
cylinder  is  at  the  origin  and  the  length  is  aligned  with  the  x-axis. 

3.1.1  Projected  Areas  from  Cylinder  Orientation 

So  let  us  apply  a  pitch-yaw-roll  rotation  sequence  to  the  cylinder  in  standard  orien¬ 
tation  and  work  out  the  projected  areas  on  the  orthogonal  y-z,  x-z,  and  x-y  planes. 
The  orientation  of  the  cylinder  is  completely  described  by  its  axis  of  rotation,  which 
is  along  the  unit  vector  k.  Let  u  be  the  final  orientation  of  the  k  vector.  Then,  using 
Table  1,  we  get 


u  =  Rfi(j)p)Rfi(j)y)k 

=  i?i(0p)(sin  0^1  -f  cos^yk) 

=  sin  (pyi  +  cos  (l>y{  —  sin  fipj  +  cos  fipk).  (3) 

It  is  also  convenient  to  specify  u  in  terms  of  the  direction  cosines: 

U  =  cos  011  -f  cos  02j  +  COS  ^sk,  (4) 


*See  Appendix  A  for  a  derivation  of  this  result. 
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where  cos  01  =  u  •  i,  cos  02  =  u  ■  j,  and  cos  03  =  u  ■  k.  Comparing  Eqs.  3  and  4 
gives  the  eorrespondenee 


01  =  cos“^(sin  (py), 

02  =  cos“^(— sin0pcos0j^),  (5) 

03  =  COS“^(cOS  0p  cos  0J,). 

Let  L  be  the  length  of  the  eylinder  and  D  be  its  diameter.  Then  the  projeeted  areas 
onto  the  orthogonal  planes  are 

TT 

LD  \  sin  01 1  +  \  cos  01 1 

LT)|  sin02|  +  COS02I .  (6) 

TT 

LD  \  sin  03 1  +  \  cos  03 1 

Therefore,  given  the  piteh  and  yaw  of  the  eylinder,  Eqs.  5  and  6  give  us  the  projeeted 
areas. 

3.1.2  Cylinder  Orientation  from  Projected  Areas 

Next,  we  work  out  the  inverse  problem  of  determining  the  piteh  and  yaw  from  the 
measured  projeeted  areas.  Thus,  given  the  length  L,  diameter  D,  and  the  projeeted 
areas  of  an  RCC  onto  the  orthogonal  planes,  we  want  to  solve  for  the  orientation  of 
the  RCC.  More  speeifieally,  we  seek  the  FATEPEN  piteh  and  yaw  rotation  sequenee 
that  will  bring  about  this  orientation.*  It  is  eonvenient  in  what  follows  to  speeify  the 
hnal  orientation  in  terms  of  the  polar  angle,  9,  measured  from  the  2:-axis  and  0,  the 
azimuthal  angle  measured  from  the  x-axis  in  the  x-y  plane: 

u  =  sin  6  cos  0i  +  sin  6  sin  0j  +  cos  Ok.  (7) 

There  is  no  loss  of  generality  if  we  restriet  u  to  lie  in  the  first  oetant  of  the  unit 
sphere,  so  that 

0  <  6  <  7r/2  and  0  <  0  <  7i/2.  (8) 


/TT 

Ayz  =  LD\  cos  1^- 

/TT 

=  LD\ cos  y- 

(TT 

A^y  =  LD\ cos 


7F  2 1  I 

n  )  I  +  -D  \  cos0i| 
+  cos  02 1 
+  jD'^l  cos  03 1 


*  Ordinarily  we  would  need  pitch,  yaw,  and  roll  to  completely  specify  orientation,  but  placing  the 
axis  of  the  RCC  along  the  z  (roll)  axis  eliminates  the  need  to  consider  roll. 
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Comparing  Eqs.  4  and  7  gives 

0  =  (j)^  and  0  =  tan“^  ( ^  _  (9) 

\COS(t)iJ 

All  the  formulas  in  Eq.  6  are  of  the  form 

asina  +  focoso;  =  Ap,  (10) 

with  a  =  (pi,  02,  or  03,  a  =  LD,  b  =  /A,  and  Ap  the  projected  area,  Ay^,  A^^, 

or  A^y.  We  can  solve  this  for  the  angle  a  by  first  introducing  another  angle  (3,  defined 
by 

0  =  tan“^(&/a),  (11) 

which  gives  a  =  y/a'^  +  b‘^  cos  0  and  b  =  y/a?  +  b‘^  sin  0,  and  then  Eq.  10  becomes 

■\/a4~+~P  cos  0  sin  a  +  y/a^  Ab"^  sin  0  cos  a  =  y/a^  Ab"^  sin(Q;  +  0)  =  Ap.  (12) 

Now  it  is  easy  to  show  that  y/a^  =  A^ax,  the  maximum  presented  area  of  the 
RCC,  so  that  0  <  Ap/ y/a^  Ab"^  <  1,  and  therefore 

sin(a  +  0)  =  .  f ^  •  (13) 

V  a  +  0 

Two  possible  solutions  to  this  equation  are 

.  -I  f  Ap  \  -1  fh\ 

a  =  sm  ,  —  tan  -  or 

\yf^?^)  \a) 

a  =  TT  —  tan“^  —  sin“^  (  ^  .  (14) 

\aj  Vv/^^TF; 

In  practice  we  try  both,  so  this  gives  2^  =  8  possible  combinations  for  0i,  02,  and 
03.  But  by  imposing  the  constraint  that 

COS^  01  +  COS^  02  +  COS^  03  =  1,  (15) 

we  will  find  that  only  one  combination  will  work.  Once  0i,  02,  and  03  have  thus 
been  found,  we  return  to  Eq.  9  to  get  9  and  0. 
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Next  we  need  a  rotation  that  will  align  k  with  u.  It  is  easy  to  check  that  Re{0),  where 


k  X  u  .  _ 

e  =  - =  —  sin  (pi  +  cos  0j, 

Ik  X  u| 


(16) 


is  this  rotation.  That  is,  Re{9)k  =  u.  The  quaternion  representation  of  this  rotation 
is 


qe{d)  =  cos(6*/2)  +  esin(6'/2)  =  cos(6*/2)  +  (—  sin  (pi  +  cos  0j)  sin(6'/2).  (17) 

Finally,  this  quaternion  can  be  factored  into  a  pitch-yaw-roll  rotation  sequence.^  The 
prescription  is  as  follows: 

1.  Set  po  =  cos(6*/2),  pi  =  —  sin(;/)sin(6'/2),  p2  =  cos  0sin(6*/2),  p^  =  0. 

2.  Set  A  =  pip2  -  PoPs,  B  =  pI- pI,  D  =  pI-  pi. 

3.  Then  =  tan“^[— 2A/(i? -f  77)]. 

4.  Set  Co  =  cos(0,./2)  and  cs  =  sin(0r/2). 

5.  Set  go  =  Po<A  -f  P3C3,  gi  =  piCo  -  P2C3,  72  =  P2C0  +  P1C3,  q^.  =  psCo  -  P0C3. 

6.  Then  =  2tan“^(gi/go)  and  0^^  =  2  tan“^(g2/go)- 

Thus,  we  accomplish  what  we  set  out  to  do:  obtain  the  pitch  and  yaw  rotation 
sequence  that  will  give  the  desired  projected  areas.  For  a  cylinder  there  are  actually 
4  orientations  that  give  the  same  projected  areas,  which  are  {(pp,(py),  {(pp,  —(py), 

(  4’p-:  4’y)-'  (  4’p-:  0J/)- 

3.2  Orientation  of  a  Cuboid 

The  cuboid,  or  RPP,  has  a  length  L,  width  W,  and  thickness  T,  where  L  >  W  >  T, 
and  is  initially  oriented  so  that  the  length  is  along  the  2:-axis,  the  width  is  along  the 
x-axis  and  the  thickness  is  along  the  y-axis.  Thus,  the  initial  areas  along  each  of  the 
axes  are  =  LT,  Ay  =  LW,  and  A^  =  WT.  We  take  i  to  be  a  unit  vector  along 
the  x-axis,  j  to  be  a  unit  vector  along  the  y-axis,  and  k  to  be  a  unit  vector  along  the 
2:-axis.  Let 

u  =  xi -I- i/j -I- ;2k  with  +  y'^  +  =  I  (18) 
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be  a  unit  vector  on  the  unit  sphere.  Then  the  projected  area  of  the  cuboid  orthogonal 
to  u  is 

Ap{x,  y,  z)  =  {AA  +  Ay]  +  A^k)  ■  u  =  A^x  +  Ayy  +  (19) 

As  we  vary  the  direction  of  u,  the  rate  of  change  of  the  projected  area  is  given  by 


DflAp  =  VAp{x,y,z)  ■  u 

=  \\V Ap{x,y,z)\\  ||u||  cos6' 

=  ||VAp(a;,|/,^)||  cos6',  (20) 


where  6  is  the  angle  between  the  gradient  and  the  unit  vector.  From  this  expression, 
it  is  clear  that  the  area  is  maximized  when  cos  6*  is  a  maximum,  which  occurs  when 
6^  =  0.  The  gradient  along  this  direction  is 

VAp{x,  y,  z)  =  =  AA  +  Ay]  +  A^k,  (21) 

and  therefore  the  magnitude  of  the  maximum  projected  area  is 

A^^=\\VAp{x,y,z)\\  =  ^Al  +  Al  +  Al,  (22) 


and  the  view  direction  for  this  maximum  is  given  by  the  unit  vector 


Ur, 


A  A  +  Ay]  + 

4  ■ 


(23) 


By  construction,  the  minimum  projected  area  is  initially  oriented  along  the  2:-axis, 
Amin  =  Az,  and  is  realized  along  the  unit  vector  Umin  =  k. 


3.2.1  Projected  Areas  from  Cuboid  Orientation 

Now  let  us  work  out  the  projected  areas  from  a  given  pitch-yaw -roll  rotation  se¬ 
quence,  R.  Again,  making  use  of  Eq.  2  and  applying  the  rotations  from  Table  1 
successively  gives 


Rj{(l)y)R^{(f)r)Axi  =  Ax  cos  0r(cos  (pyi  —  sin  ipyk)  +  Ax  sin  (pA, 
RA(l>y)Ri.{(l>r)Ay]  =  —Ay  sin  ^^(cos  (pyi  —  sin  (pyk)  +  Ay  cos  <pr],  (24) 

RA<Py)Ri.{<Pr)Azk  =  Azismcpyl  +  cos^yk). 
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Collecting  terms  we  get 


Rl{(l)y)R^{(j)r)Axi  =  cos  (pr  cos  (pyi  +  Ay;  sin  pri  —  A^  cos  pr  sin  pyk, 
Rl{(l)y)R^{(l)r)Ay}  =  —Ay  sm  (pr  cos  pyi  +  Ay  COS  pr3  +  Ay  sin  pr  sin  pyk,  (25) 
R~^{py)R^{pr)Azk  =  Az  sin py\  +  A^  cos  pyk. 

Finally,  applying  pitch,  we  get 

RAA  =  A  Ay.  cos  pr  cos  pyi 

+  Ay  sin  pr{cos  0pj  +  sin  0pk) 

—  Ay  COS  pr  sin  py{—  sin  +  cos  ppk), 

RAy]  =  —  Ay  sin  pr  cos  py\  (26) 

+  Ay  COS  pricos  ppl  +  sin  pp}^ 

+  Ay  sin  pr  sin  py{—  sin  0pj  +  cos  0pk), 

RAyk  =  +  Ay  sin  pyi 

+  Ay  COS  py{—  sin  0pj  +  cos  ppk). 


and  collecting  terms, 

RAyi  =  +  Ay  cos  pr  cos  pyi 

+  Aj.(sin  pr  cos  pp  +  cos  pr  sin  py  sin  pp)} 

+  A2;(sin  pr  sin  pp  —  cos  pr  sin  py  cos  0p)k,  (27) 

RAyl  =  —  Ay  sin  pr  cos  pyi 

+  Ap(cos  pr  COS  Pp  —  sin  pr  sin  py  sin  pp)] 

+  ylp(cos  pr  sin  Pp  +  sin  pr  sin  py  cos  pp)k,  (28) 

RAyk  =  +  Ay  sin  pyi  —  Ay  cos  py  sin  0pj  +  Az  cos  py  cos  ppk.  (29) 

Now  let  Ayy  be  the  projection  on  the  x-y  plane,  Ayz  be  the  projection  on  the  x-z 
plane,  and  Ayz  be  the  projection  on  the  y-z  plane.  Then  we  have 

Ayy  =  [R{Ayi  +  Ay}  +  ^zk)]  •  k, 

Ayz  =  [R{Ayi  +  Ay}  +  A^k)]  ■  j, 

Ayz  =  [R{,Ayi  +  Ay}  +  A^k)]  ■  i. 
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(30) 

(31) 

(32) 


Three  is  the  maximum  number  of  sides  that  can  project  unto  a  plane  for  a  cuboid,  so 
we  can  include  all  the  contributions  by  using  the  absolute  value.  Thus,  we  get  the 
projected  areas 

Axy  =Ax  I  sin  (pr  sin  —  cos  (pr  sin  <py  cos  (pp  \  + 

Ay  I  cos  (pr  sin  (pp  +  sin  (pr  sin  (py  cos  (pp  \  + 

Az  \  cos  <py  cos  (pp\,  (33) 

A^z  =Ax  I  sin  (pr  cos  (pp  +  cos  (pr  sin  (py  sin  (pp  \  + 

Ay  I  COS  (pr  COS  (pp  —  sin  (pr  sin  (py  sin  (pp  \  + 

Az  I  cos  <py  sin  (pp  I ,  (34) 

Ayz  =Ax  I  cos  (pr  cos  (py  \  +  Ay  \  sin  (pr  cos  (py  \  +  Az  \  sin  (py\.  (35) 


3.2.2  Cuboid  Orientation  from  Projected  Areas 

Next,  given  measured  values  for  A^y,  A^z,  and  Ayz,  we  want  to  solve  Eqs.  33-35 
for  (pp,  (py,  and  (pr.  This  is  a  system  of  3  simultaneous  nonlinear  equations,  which 
can  be  solved  numerically  by  Newton’s  method  in  3  dimensions.  Define 


f{(pp,  (py,  (pr)  =Ax\  sinc))^  sin0p  —  cos(pr  sin^j^  cos0p|  + 

Ay  I  cos  (pr  sin  (pp  +  sin  <pr  sin  (py  cos  (pp  \  + 

cos  (py  cos  (pp\  -  Axy,  (36) 


g{(pp,  (py,  (pr)  =Ax\  sin0r  cos  (pp  +  cos  (pr  sin  (py  sin  (pp\  + 

Ay  I  cos  (pr  cos  (pp  —  sin  (pr  sin  (py  sin  (pp  \  + 

A2I  cos  (py  sin  (pp\  -  Axz,  (37) 

h{(pp,  (py,  (pr)  =  Ax  \  cos  (pr  COS  (py  \  +  Ayj  sin^r  cos  (py  \  +  Az  \  sin  02^ I  —  Ayz.  (38) 


Then  the  problem  is  to  find  the  roots  of  /,  g,  and  h.  Let  A  be  the  matrix  of  partial 
derivatives 


f<j>p{.4’p)4’yi4’r)  f 4’yi.4’p-i  4’y)  4’r)  f 4>ri.4’p)  4’y)  4’r) 

9<i>p{4^pi 4^y> 4^r)  g4>yi.4^pi 4^yi 4^r)  gcpri^^pi 4^r) 

}^<l>p{.4’pi4’yi4’r)  h(j)y{(pp,  (py,  (pr)  h^^{(pp,  (py,  (pr) 


(39) 
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where 


(40) 


f(t>p{4^pi  4^y>  4^r) 


\4^p')  4^y)  4^r 

d<f)p 


and  so  on  for  the  other  partial  derivatives.  For  this  we  need  the  derivative  of  the 
absolute  value: 


uu'  uu' 


d  ,  ,  d  !  ^  1  .  o.  1  /  Li  Li  Li  Li  /  \  /  .  ,  ^ 

—  \u\  =  —vv?  =  -iu^)  ^^2uu  =  =  — —  =  Sgn(M)  M  ,  (41) 

dx  dx  2  ^/u^  \u\ 

where  sgn(a;)  is  the  sign  function: 


sgn(a;)  =  < 


+1  if  a;  >  0 
0  if  a;  =  0  , 
—  1  if  a;  <  0 


(42) 


Then  from  Eqs.  36-38  we  get 


f<l>p  i.4’p)  4’y)4’r)  =  +^x  sgn(sin  (pr  sin  (pp  —  cos  (pr  sin  (py  cos  ( 

/  •  /  /  I  I  •  I  •  I  \ 


f4’yi.4’pi4’yi4’r) 


frt>r{^Pi4^yi4^r) 


(sin  (pr  cos  (pp  +  cos  (pr  sin  (py  sin  ippj 
+Ay  sgn(cos  (pr  sin  (pp  +  sin  (pr  sin  (py  cos  ( 
(cos  (pr  cos  (pp  —  sin  (pr  sin  (py  sin  (pp) 
+Az  sgn{cos  (py  cos  (pp){—  cos  (py  sin0p), 

+Ax  sgn(sin  (pr  sin  (pp  —  cos  (pr  sin  (py  cos  ( 
(—  cos  (pr  cos  (py  sin  (pp) 

+Ay  sgn(cos  (pr  sin  (pp  +  sin  (pr  sin  (py  cos  ( 
(sin  (pr  cos  (py  cos  (pp) 

+Az  sgn(cos0j^  cos  (pp){—  sin  (py  cos(pp), 

+Ax  sgn(sin  (pr  sin  (pp  —  cos  (pr  sin  (py  cos  ( 
(cos  (pr  sin  (pp  +  sin  (pr  sin  (py  cos  (pp) 
+Ay  sgn(cos  (pr  sin  (pp  +  sin  (pr  sin  (py  cos  ( 
(—  sin  (pr  sin  (pp  +  cos  (pr  sin  (py  cos  0p), 


(43) 


(44) 


(45) 
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9(l>p{^pi  4^yi  4^r) 


9<t>y  4^yi  4^r) 


-  COS  (pr  sin  (pp 
+Az  sgn(cos  (py  sin  (pp)  (cos  (py  cos  <pp) , 


9(f>r  {4^pi  4^y^ 


+Ax  sgn(sin  (pr  cos  (pp  +  cos  (pr  sin  (py  sin  cp. 

(—  sin  (pr  sin  (pp  +  cos  (pr  sin  (py  cos  (pp) 
+Ay  sgn(cos  (pr  cos  (pp  —  sin  (pr  sin  (py  sin  0, 

(py  cos  (pp\ 

'cos 

+Ax  sgn(sin 
(cos  (pr  COS  (py  sin  i^p 
+Ay  sgn(cos  (pr  cos  (pp  —  sin  (pr 

r  COS  (py  sin  (pp 

+A^  sgn(( 

+Ay;  sgn(sin  (pr  cos  (pp  +  cos  (pr  sin  (py  sin 
(pr  sin  (py  sin  (pp 


sin  (pr  sin  ( 

(cos  (py 

.  (pr  COS  (pp  +  COS  (pr  sin  (py  sin  ( 
th..  ‘^in  thA 


sin  (py  sin  (f.pj 

sin  (pr  cos  (py  sin  ( 

cos  (py  sin  (pp){—  sin  (py  sin  (pp), 

sin  (py  s 


(cos  (pr  cos  (pp  —  sin  , 

+Ay  Sgn(cOS  (pr  cos  (pp 
i —  ciri  /A  /A  —  nr 


sin  (pr  sin  (py 

—  sin  (pr  cos  (pp  —  cos  (pr  sin  (py 


sine 
sin0p), 


hrj>y  {(ppi  (py,  (pr 


)=  0, 


:  +Ax  Sgn(cOS  (pr  cos  (py){—  cos  (pr  sin0y)  + 
+Ay  sgn(sin  (pr  cos  (py){—  sin  (pr  sin0y)  + 
+A,,  sgn(sin  epy) (cos  (py), 

Kr  {(pp^  (py:  (pr)  =  +(^x  Sgn(cOS  (pr  cos  (py){- sin  (pr  cos  (py) 
+Ay  sgn(sin  (pr  cos  (py) (cos  (pr  cos  (py). 


(46) 


(47) 


(48) 

(49) 


(50) 


(51) 


Newton’s  method  for  solving  this  3-dimensional  (3D)  problem  numerically  for  pitch, 
yaw,  and  roll  is  expressed  by  the  matrix  iteration  equation 


(52) 


where  the  matrix  A  is  given  by  Eq.  39.  We  pick  a  starting  orientation  and  continue 
iterating  until  either  the  (5’s  fall  below  a  preset  tolerance  value  e  or  we  exceed  a 
maximum  number  of  iterations — in  which  case  we  try  another  starting  orientation. 


d(pp 

0p,n+l 

f  {(pp:  (py:  (pr) 

S(Py 

= 

4^y,n 

=  A-^ 

9{(pp:  (py:  (pr) 

6(pr 

_0r,n+l 

4^r,n_ 

h{(pp:  (py,  (pr) _ 
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23 

24 
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26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 
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41 

42 
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45 

46 

47 

48 

49 
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If  we  write  this  eoeffieient  matrix  (Eq.  39)  as 


Oil 

0-12 

Ol3 

A  = 

<321 

0-22 

<323 

<331 

O32 

<333 

then  the  inverse  is  given  by 


(53) 


1 

det  A 


O22O33  ~  023^32 
^23031  ~  0,210^33 
<^21(^32  —  ^22031 


013^32  —  <312033  ai2<323 

<311033  ~  013031  013021 

O12O31  —  O11O32  O11O22 


O13O22 

O11O23 

O12O2I 


(54) 


Thus,  we  now  have  all  the  ingredients  for  implementing  the  proposed  solution  with 
the  program  in  Listing  1 . 

Listing  1.  orient. cpp 

//  orient. cpp:  Given  the  dimensions  of  an  RPP  and  projected  areas  onto  three  orthogonal  planes, 


// 

finds  a  pitch -yaw- 

roll 

rotation 

sequence  that 

will 

orient  the  RPP  for  these  projected  areas. 

1/ 

The  orientation  is 

not 

unique  and  will  depend 

upon 

the  initial  starting  point. 

// 

Method  of  solution 

makes  use  of 

Newton's  method  in 

three  dimensions. 

//  R.  Saucier,  October  2015 

#include  <iostream> 

#include  <cstdlib> 

#include  <cmath> 

#include  <iomanip> 

#include  <chrono> 

#include  <random> 

inline  double  sgn(  double  x  )  { 

if  {  X  >  0.  )  return  +1. ; 
else  if  {  X  <  0.  )  return  -1.; 
else  return  0.; 


int  main(  int  argc,  char*  argv[]  )  { 

const  double  D2R  =  M_PI  /  180.;  //  to  convert  deg  to  rad 

const  double  R2D  =  180.  /  PEPI;  //  to  convert  rad  to  deg 

const  double  TOL  =  l.e-9;  //  convergence  criterion 

const  int  N  =  100;  //  max  number  of  iterations 


double 

L  =  3. ; 

// 

length 

double 

W  =  2. ; 

// 

width 

double 

T  =  1. ; 

// 

thickness 

double  Ax  =  L  *  T;  //  area  of  side  (initially  along  i) 

double  Ay  =  L  *  W;  //  area  of  top  (initially  along  j) 

double  Az  =  W  *  T;  //  area  of  front  (initially  along  k) 

double  A_xy  =  4.5;  //  projected  area  on  x-y  plane 

double  A_x2  =  6.5;  //  projected  area  on  x-z  plane 

double  A_yz  =  3.5;  //  projected  area  on  y-z  plane 

double  p,  y,  r,  cp,  sp,  cy,  sy,  cr,  sr,  f,  g,  h; 

double  el,  e2,  e3,  e4,  e5,  e6,  e7,  e8,  e9; 

double  si,  s2,  s3,  s4,  s5,  s6,  s7,  s8,  s9; 

double  fp,  fy,  fr,  gp,  gy,  gr,  hp,  hy,  hr; 

double  all,  al2,  al3,  a21,  a22,  a23 ,  a31,  a32,  a33,  det; 

double  bll,  bl2,  bl3,  b21,  b22,  b23 ,  b31,  b32,  b33; 

double  del_p,  del_y,  del_r; 

//  initial  estimate  can't  be  zero 

unsigned  int  seed  =  std :: chrono :: high_resolution_clock :: now( ). time_since_epoch (). count () ; 
std::mtl9937  rng(  seed  );  //  Mersenne  Twister  engine 
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51 

std : : uniform- real- 

distribution<double>  uniform(  0.,  M-PI  );  //  uniform  distribution 

52 

p  =  uniform(  rng 

; 

53 

y  =  uniform(  rng 

; 

54 

r  =  uniform(  rng 

; 

55 

56 

if  (  argc  ==  4  )  {  //  optionally  specify  initial  pitch,  yaw,  roll  on  commandline 

57 

58 

P  = 

atof(  argv[l]  )  *  D2R; 

59 

y  = 

atof(  argv[2]  )  *  D2R; 

60 

r  = 

atof(  argv[3]  )  *  D2R; 

61 

> 

62 

63 

for  ( 

int  i  =  0;  i  <  N;  i++  )  { 

64 

65 

cp 

=  cos(  p  ) ; 

sp  =  sln(  p  ) ; 

66 

cy 

=  cos(  y  ) ; 

sy  =  sln(  y  ) ; 

67 

cr 

=  cos(  r  ) ; 

sr  =  sin (  r  ) ; 

68 

69 

el 

=  sr=t=sp- 

cr  *  sy  *  cp; 

70 

e2 

=cr*sp+ 

sr  *  sy  *  cp; 

71 

e3 

=  cy  *  cp; 

72 

e4 

=sr*cp+ 

cr  *  sy  *  sp; 

73 

e5 

=  cr=t=cp- 

sr  *  sy  *  sp; 

74 

e6 

=  cy  *  sp; 

75 

e7 

=  cr  *  cy; 

76 

e8 

=  sr  *  cy; 

77 

e9 

=  sy; 

78 

79 

si 

=  sgn(  el  ) 

80 

s2 

=  sgn(  e2  ) 

81 

S3 

=  sgn(  e3  ) 

82 

s4 

=  sgn  e4  ) 

83 

s5 

=  sgn(  e5  ) 

84 

s6 

=  sgn(  e6  ) 

85 

s7 

=  sgn(  e7  ) 

86 

s8 

=  sgn(  e8  ) 

87 

s9 

=  sgn(  e9  ) 

88 

89 

f  = 

Ax  *  fabs( 

el  )  +  Ay  *  fabs(  e2  )  +  Az  *  fabs(  e3  )  -  A-xy; 

90 

g  = 

Ax  *  fabs( 

e4  )  +  Ay  *  fabs(  e5  )  +  Az  *  fabs(  e6  )  -  A-xz; 

91 

h  = 

Ax  *  fabs( 

e7  )  +  Ay  *  fabs(  e8  )  +  Az  *  fabs(  e9  )  -  A-yz; 

92 

93 

fp 

=  Ax  =t=  si  * 

(  sr  *  cp  +  cr  *  sy  *  sp  )  +  Ay  *  s2  *  (  cr  *  cp  -  sr  *  sy  *  sp  )  +  Az  * 

s3  *  (  -cy  *  sp  ) ; 

94 

fy 

=  Ax  *  si  * 

(  -cr  *  cy  ♦  sp  )  +  Ay  *  s2  *  {  sr  *  cy  *  cp  )  +  Az  *  s3  *  (  -sy  *  cp  ) ; 

95 

fr 

=  Ax  =t=  si  * 

(  cr  *  sp  +  sr  *  sy  *  cp  )  -1-  Ay  *  s2  *  (  -sr  *  sp  +  cr  *  sy  *  cp  ) ; 

96 

97 

gp 

=  Ax  *  s4  * 

(  -sr  *  sp  -1-  cr  *  sy  *  cp  )  +  Ay  *  s5  *  (  -cr  *  sp  -  sr  *  sy  *  cp  )  +  Az 

*  s6  *  (  cy  *  cp  ) ; 

98 

gy 

* 

X 

< 

(  cr  *  cy  *  sp  )  -1-  Ay  *  s5  *  (  -sr  *  cy  *  sp  )  -i-  Az  *  s6  *  (  -sy  *  sp  ) ; 

99 

gr 

=  Ax  *  s4  * 

(  cr  *  cp  -  sr  *  sy  ♦  sp  )  -1-  Ay  *  s5  *  (  -sr  *  cp  -  cr  *  sy  *  sp  ) ; 

100 

101 

hp 

=  0.; 

102 

hy 

=  Ax  *  s7  * 

(  -cr  *  sy  )  -1-  Ay  *  s8  *  (  -sr  *  sy  )  -i-  Az  *  s9  *  {  cy  ) ; 

103 

hr 

=  Ax  *  s7  * 

(  -sr  *  cy  )  -1-  Ay  *  s8  *  (  cr  *  cy  ) ; 

m 

105 

all 

=  fp;  al2 

=  fy;  al3  =  fr; 

106 

a21 

=  gp;  a22  - 

=  gy;  a23  =  gr; 

107 

a31 

=  hp;  a32 

=  hy;  a33  =  hr; 

108 

det 

=  all  *  (  a22  *  a33  -  a23  *  a32  )  +  al2  *  (  a23  *  a31  -  a21  *  a33  )  +  al3  *  (  a21  * 

a32  -  a22  *  a31  ) ; 

109 

if 

(  det  ==  0 

{ 

no 

111  std::cerr  «  "bad  initial  orientation:  "  «  p  *  R2D  «  "\t"  «  y  *  R2D  «  "\t"  «  r  *  R2D  «  std::endl 

112  «  "program  "  «  argv[0]  «  "  stopped"  «  std::endl; 

113  exit(  EXIT-FAILURE  ); 

114  } 

115 


116 

bll  = 

(  a22 

* 

a33 

-  a23 

♦ 

a32  ) 

/ 

det ; 

117 

bl2  = 

(  al3 

* 

a32 

-  al2 

* 

a33  ) 

/ 

det ; 

118 

bl3  = 

(  al2 

* 

a23 

-  al3 

* 

a22  ) 

/ 

det ; 

119 

b21  = 

(  a23 

* 

a31 

-  a21 

* 

a33  ) 

/ 

det ; 

120 

b22  = 

(  all 

* 

a33 

-  al3 

* 

a31  ) 

/ 

det ; 

121 

b23  = 

(  al3 

* 

a21 

-  all 

* 

a23  ) 

/ 

det ; 

122 

b31  = 

(  a21 

* 

a32 

-  a22 

* 

a31  ) 

/ 

det ; 

123 

b32  = 

(  al2 

* 

a31 

-  all 

* 

a32  ) 

/ 

det ; 

124 

b33  = 

(  all 

* 

a22 

-  al2 

* 

a21  ) 

/ 

det ; 

125 

126 

del-p 

=  bll 

* 

f  + 

bl2  * 

g 

-1-  bl3 

* 

h ; 

127 

del-y 

=  b21 

* 

f  + 

b22  * 

g 

-1-  b23 

* 

h; 

128 

del-T 

=  b31 

* 

f  + 

b32  * 

g 

-1-  b33 

* 

h ; 

129 

130  p  -=  del_p: 

131  y  -=  del_y: 

132  r  -=  del_r: 

133 

134  if  (  fabs(  del_p  )  <  TOL  &&  fabs(  deUy  )  <  TOL  &&  fabs(  del_r  )  <  TOL  )  break; 

135 
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136 

137 

138 

139 

140 

141 

142 

143 

144 

145 

146 

147 

148 

149 

150 

1 

2 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 


if  (  i  ==  N-1  )  { 


std::cerr  «  "failed  to  converge  after  "  «  N  «  "  iterations:  f  =  "  «  f  «  "  g  =  "  «  g  «  "  h  =  "  «  h  «  std 
: ; endl 

«  "try  another  initial  orientation"  «  std::endl 
«  "program  "  «  argv[0]  «  "  stopped"  «  std::endl; 
exit(  EXIT-FAILURE  ); 

} 

> 

std : : cout  «  std : : setprecision (9)  «  std::fixed; 
std : : cout  «  "pitch  =  "  «  fmod(  p  *  R2D,  360.  )  «  std:: endl 

«  "yaw  =  "  «  fmod(  y  *  R2D,  360.  )  «  std:: endl 

«  "roll  =  "  «  fmod(  r  *  R2D,  360.  )  «  std:: endl; 

return  EXIT-SUCCESS; 


The  program  may  be  eompiled  and  run  with  the  following  eommands: 


C++  -02  -Wall  -std=c++ll  -o  orient  orient. cpp  -Im 
./orient 


4.  Results  and  Discussion 

We  now  have 

•  a  method  for  computing  the  projected  areas  on  orthogonal  planes,  given  the 
cylinder  or  cuboid  orientation,  and 

•  a  method  for  determining  the  pitch-yaw-roll  rotation  sequence  to  bring  about 
the  orientation,  given  the  projected  area  measurements. 

This  means  that  we  have  a  way  to  check  the  results.  We  will  show  that  projected  area 
measurements  lead  to  a  unique  orientation  for  a  cylinder,  but  that  is  not  the  case  for 
a  cuboid. 

4.1  Cylinder  Orientation  Sample  Case 

Consider  an  RCC  with  L  =  1,  D  =  1,  (pp  =  30°,  and  (py  =  15°.  The  projected  areas 
are  found  to  be  Ay^  =  1.16920,  =  1.25496,  and  A^y  =  1.20494.  Using  these 

values  for  the  projected  areas,  the  program  in  Listing  2  computes  (pp  =  ±30°  and 
(py  =  ±15°,  with  the  4  possible  orientations  displayed  in  Fig.  2. 

Listing  2.  rcc.cpp 

//  rcc.cpp;  compute  pitch  and  yaw  from  area  projections  for  an  RCC  (sample  case) 

//  R.  Saucier,  June  2006  (Revised  October  2015) 

#include  "Rotation. h" 

#include  <iostream> 

#include  <cmath> 

#include  <cstdlib> 

#include  <cassert> 

#include  <iomanip> 
using  namespace  std; 

double  anglel{  double  a,  double  b,  double  c  )  {  return  asin(  c  /  sqrt(  a  *  a  +  b  +  b  )  )  -  atan(  b  /  a  );  } 
double  angle2(  double  a,  double  b,  double  c  )  {  return  M-PI  -  asin(  c  /  sqrt(  a  *  a  +  b  +  b  )  )  -  atan{  b  /  a  );  } 

int  main(  void  )  { 
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16 

17 

18 

19 

20 
21 
22 

23 

24 

25 


32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 
61 
62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 


const  double  L  =  1.,  0=1.; 

const  double  A  =  L  *  D; 

const  double  B  =  0.25  *  I^PI  *  D  *  D; 

const  double  A_MIN  =  A  <  B  ?  A  ;  B; 

const  double  A_MAX  =  sqrt(  A  *  A  +  B  *  B 

cout  «  setprecision (6)  «  fixed; 


cout  «  "L 
cout  «  "D 


«  L  «  endl ; 
«  D  «  endl ; 


26 

cout  «  "A_MIN  =  "  « 

A_MIN 

«  endl : 

27 

cout  «  "A_MAX  =  "  « 

A_MAX 

«  endl : 

28 

29 

double  A_xy  =  1.20494 

// 

projected 

area 

on 

x-y 

plane 

30 

double  A_xz  =  1.25496 

// 

projected 

area 

on 

x-z 

plane 

31 

double  A_yz  =  1.16920 

// 

projected 

area 

on 

y-z 

plane 

cout  «  "A_xy  = 
cout  «  "A_X2  = 
cout  «  "A_yz  = 


«  A_xy  «  endl; 
«  A_xz  «  endl; 
«  A_yz  «  endl; 


assert (  A_MIN  <=  A_xy  &&  A_xy  <=  A_MAX  ) 

assert (  A_MIN  <=  A_yz  &&  A_yz  <=  A_MAX  ) 

assert (  A_MIN  <=  A_xz  &&  A_xz  <=  A_MAX  ) 

double  phi_x[  2  ],  phi_y[  2  ],  phi_z[  2  ]; 
phi_x[0]  =  anglel{  A,  B,  A_yz  ); 

phi_x[l]  =  angle2(  A,  B,  A_yz  ); 

phi_y[0]  =  anglel(  A,  B,  A_xz  ); 

phi_y[l]  =  angle2{  A,  B,  A_xz  ); 

phi_z[0]  =  anglel{  A,  B,  A_xy  ); 

phi_z[l]  =  angle2(  A,  B,  A_xy  ); 

double  d,  delta  =  l.e36; 
int  ii  =  0,  jj  =0,  kk  =  0; 
for  (  int  i  =  0;  i  <  2;  i++  )  { 
for  (  int  j  =0;  j  <2;  j++  )  { 
for  (  int  k  =  0;  k  <  2  ;  k++  )  { 

d  =  cos(  phi_x[i]  )  *  cos(  phi_x[i]  )  + 

cos(  phi_y[j]  )  *  cos(  phi_y[j]  )  + 

cos(  phi_z[k]  )  *  cos(  phi_z[k]  ); 

if  (  fabs(  d  -  1.  )  <  delta  )  { 

ii  =  i; 

jj  =  j; 

kk  =  k; 

delta  =  fabs(  d  -  1.  ) ; 


double  th  =  phi_z[kk]; 

double  ph  =  atan(  cos(  phi_y[jj]  )  /  cos(  phi_x[ii] 


73 

cout  « 

"Derived  Angles 

(deg):  " 

«  endl 

74 

« 

"  phi_x: 

"  « 

phi_x[ii; 

*  va : : R2D  « 

endl 

75 

« 

phi_y : 

"  « 

Phi-y[jj] 

*  va : : R2D  « 

endl 

76 

« 

"  phi_2: 

"  « 

phi_z[kk; 

*  va : : R2D  « 

endl 

77 

« 

"  theta: 

"  « 

th  ♦  va : : 

R2D  «  endl 

78 

« 

"  phi : 

"  « 

ph  *  va : : 

R2D  «  endl; 

80 

81 

82 

83 

84 

85 

86 

87 

88 

89 

90 

91 

92 

93 

94 

95 


va::Vector  u  =  sin(  th  )  *  cos{  ph  )  *  va::Vector(  1.,  0.,  0.  )  + 
sin(  th  )  *  sin{  ph  )  *  va::Vector(  0.,  1.,  0.  )  + 
cos(  th  )  *  va::Vector(  0.,  0.,  1.  ); 

va::Rotation  R(  va::\/ector(  0. ,  0. ,  1.  ) ,  u  ) ;  //  find  rotation  that  takes  k  to  u 

va::sequence  s  =  factor(  R,  va::XYZ  ); 


cout 

« 

"pyr  = 

« 

+s . first 

* 

va : 

:R2D 

« 

"Xt" 

« 

+s . second 

va : 

:R2D 

« 

,.\t" 

« 

s. third 

va : 

:R2D 

« 

endl ; 

cout 

« 

‘•pyr  = 

« 

+s . first 

* 

va : 

:R2D 

« 

"\t" 

« 

-s . second 

va : 

:R2D 

« 

"Xt" 

« 

s . third 

va : 

:R2D 

« 

endl ; 

cout 

« 

"pyr  = 

« 

-s. first 

* 

va : 

:R2D 

« 

"\t" 

« 

+s . second 

va : 

:R2D 

« 

"Xt" 

« 

s. third 

va : 

:R2D 

« 

endl ; 

cout 

« 

■•pyr  = 

« 

-s. first 

* 

va : 

.•R2D 

« 

"\t" 

« 

-s . second 

va : 

:R2D 

« 

"Xt" 

« 

s . third 

va : 

:R2D 

« 

endl ; 

cout  «  "Total  effective  yaw  =  "  «  acos(  cos(  s. first  )  *  cos(  s. second  )  )  *  va::R2D  «  endl; 


return  EXIT_SUCCESS; 

} 
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Fig.  2.  An  oriented  cylinder  for  (j)p  =  ±30°  and  (/)j,  =  ±15°  is  shown  with  its  projected  area  in  the 
x-y  plane  in  dark  blue,  the  x-z  plane  in  green,  and  the  y-z  plane  in  red.  All  of  these  have  the 
same  projected  areas  on  each  of  the  orthogonal  planes.  The  actual  orientations  are  not  quite 
the  same,  but  these  are  superficial  differences  since  the  effective  yaw  angle  and  the  impact 
angle  as  dehned  by  FATEPEN  are  all  the  same  (33.2259°  in  this  example)  and  thus  give  the 
same  results  for  plate  impact. 

4.2  Cuboid  Orientation  Sampie  Case 

Consider  an  RPP  with  L  =  3  cm,  W  =  2  cm,  and  T  =  1  cm,  and  let  the  projected 
areas  be  A^y  =  4.5  cm^,  A^z  =  6.5  cm^,  Ay^  =  3.5  cm^.  When  this  information  is 
given  to  the  orient .  cpp  program  (Listing  1),  it  hnds  4  different  orientations  that  will 
work,  as  displayed  in  Fig.  3,  demonstrating  that  the  orientation  is  not  unique. 


Eig.  3.  The  projected  area  on  the  orthogonal  planes  does  not  uniquely  determine  the  cuboid  orienta¬ 
tion.  The  oriented  cuboid  (L  =  3  cm,  W  =  2  cm,  T  =  1  cm)  is  shown  with  its  projected  area 
in  the  x-y  plane  in  dark  blue  (4.5  cm^),  the  x-z  plane  in  green  (6.5  cm^),  and  the  y-z  plane  in 
red  (3.5  cm^).  All  4  have  the  same  projected  areas  on  each  of  the  orthogonal  planes,  although 
we  see  that  the  actual  orientations  are  not  the  same.  One  consequence  of  this  is  that  the  impact 
angle,  fi,  which  is  the  minimum  angle  that  a  face  makes  with  the  x-y  plane,  along  with  the 
effective  yaw  angle,  have  different  values  for  the  4  orientations — and  this  will  lead  to 
different  predictions  from  EATEPEN. 


Next,  we  want  to  see  what  difference,  if  any,  these  different  orientations  will  make 
in  FATEPEN.  The  cuboid  penetrator  is  taken  to  be  steel  4140  with  a  Brinell  hardness 
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of  300.  It  has  a  length  of  3  cm  =  1.1811  inches,  width  of  2  cm  =  0.7874  inch,  and 
thickness  of  1  cm  =  0.3937  inch.  Using  a  steel  density  of  7.83  g/cm^,  it  has  a  mass  of 
725  gr.  The  plate  target  is  taken  to  be  0.25-inch  steel  with  a  Brinell  hardness  of  150. 
All  impacts  were  at  0°  obliquity.  The  cuboid  is  initially  aligned  so  that  its  length 
is  along  the  2;-axis,  its  width  along  the  x-axis,  and  it  thickness  along  the  y-axis. 
Then  it  is  given  the  pitch-yaw-roll  rotation  sequence  listed  in  Table  2  to  give  it  the 
orientation  shown  in  Fig.  3,  resulting  in  the  following  projected  areas:  A^y  =  4.5 
cm^,  =  6.5  cm^,  and  Ay^  =  3.5  cm^.  The  target  is  taken  to  lie  in  the  x-y  plane. 
This  table  also  lists  the  residual  mass  rrir,  residual  velocity  Vr,  and  limit  velocity  vl 
for  the  4  orientations.  We  see  that  there  are  significant  differences  in  the  residual 
velocities.  It  is  important  to  emphasize  that  all  4  orientations  have  the  same  projected 
areas. 


Table  2.  FATEPEN  predictions  for  a  725-gr  steel  cuboid  with  a  striking  velocity  of  1500  f/s  impacting 
a  0.25-inch  mild  steel  plate  using  version  3.3.10.3.  In  each  case,  the  projected  areas  on  the 
orthogonal  planes  are  the  same.  Note,  in  particular,  the  sensitivity  of  the  residual  velocity,  Vr- 


Pitch  (deg) 

Yaw  (deg) 

Roll  (deg) 

fi  (deg) 

(deg) 

nir  (gr) 

Vr  (f/s) 

Vl  (f/s) 

-160.576322 

193.931472 

1.054946 

23.7 

23.7 

720.341 

841.729 

1194.12 

-166.139458 

299.729320 

5.443168 

61.2 

61.2 

723.672 

547.712 

1345.23 

176.987194 

67.658507 

13.759713 

67.7 

67.7 

723.349 

376.671 

1424.57 

205.968278 

176.788266 

183.822078 

26.2 

26.2 

721.834 

830.670 

1195.00 

4.2.1  Impact  Angle  and  Effective  Yaw  Angle 

Let  us  also  document  how  the  impact  angle  and  effective  yaw  angles  are  calculated 
to  check  on  the  values  output  by  FATEPEN.  Eet  us  consider  the  area  vector  for  the 
front,  side,  and  top  of  the  cuboid.  Initially, 

Af  =  WTk,  A,  =  LTi,  and  A  =  LWl  (55) 

And  after  the  pitch-yaw-roll  rotation  sequence,  these  vectors  can  be  calculated  by 
adapting  Eqs.  30-32: 

A'^  =  1UT[  sin  0^,1  —  cos  sin  0pj  -F  cos(f)y  cos^pk],  (56) 

A'  =  LT[  (cos^r  cos0p)i-|- 

(sin  (pr  cos  (pp  +  cos  (pr  sin  (py  sin  (pp)]+ 

(sin  (pr  sin  (pp  —  cos  (pr  sin  (py  cos  0p)k],  (57) 
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(58) 


=  LW[{—  sin  (pr  cos  0j,)i+ 

(cos  (pr  cos  pp  —  sin  pr  sin  py  sin  0p)j+ 

(cos  pr  sin  pp  +  sin  pj.  sin  py  cos  0p)k] . 

The  effective  yaw  angle,  p^s,  is  the  angle  between  the  length  of  the  euboid  and  the 
target  normal.  It  is  obtained  from  the  dot  produet  between  and  k,  so  from  Eq.  56, 
0eff  =  cos“^  (cos  py  cos  pp) .  However,  if  this  angle  turns  out  to  be  greater  than  90°, 
then  we  should  use  the  baek  faee.  Therefore,  the  eorreet  angle  is  aetually 

^eff  =  min[cos“^(cos(;/)j^cos(;/)p),  tt  —  cos“^(cos0y  cos^p)].  (59) 

The  contact  angle,  /i,  is  the  angle  between  eaeh  faee  of  the  euboid  and  the  target 
normal,  and  the  impact  angle,  p,  is  the  minimum  of  all  the  eontaet  angles.  Therefore, 
from  Eqs.  55,  56,  and  57, 

p  =  min[cos“^(coS(;/)p  cos^p), 

TT  —  COS~^  {cospy  COSpp), 

cos“^(sin  pr  sinpp  —  cos  pr  sin  py  cospp), 

TT  —  cos“^(sin  pr  sinpp  —  cos  pr  sin  py  cospp), 
cos~^ {cos pr  sinpp  +  sinpr  sin  py  cospp), 

TT  —  cos“^(cos(;/)r  sin^p  +  sin^,,  sin  py  cos0p)].  (60) 

Applying  these  formulas  for  the  4  orientations  shown  in  Table  2,  we  find  p  = 
23.7454°,  p  =  30.4564°,  p  =  24.5429°,  and  p  =  26.1524°.  Notiee  that  2  of  these 
angles  disagree  with  the  values  output  by  EATEPEN  as  listed  in  Table  2.  Even  though 
the  EATEPEN  doeumentation  states  that  p  is  the  minimum  eontaet  angle,  it  seems  to 
always  equal  the  effeetive  yaw  angle  instead. 

A  newer  version  of  EATEPEN  (3.3.18.0)  was  also  run,*  whieh  produeed  the  results 
listed  in  Table  3. 


*I  thank  Timothy  Mallory  for  making  rnns  with  both  version  3.3.16.10  and  3.3.18.0,  which 
produced  the  same  results  shown  in  Table  3. 
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Table  3.  FATEPEN  predictions  for  a  725-gr  steel  cuboid  with  a  striking  velocity  of  1500  f/s  impacting 
a  0.25-inch  mild  steel  plate  using  version  3.3.18.0.  In  each  case,  the  projected  areas  on  the 
orthogonal  planes  are  the  same. 


Pitch  (deg) 

Yaw  (deg) 

RoU  (deg) 

fi  (deg) 

(deg) 

rur  (gr) 

Vr  (f/s) 

Vl  (f/s) 

-160.576322 

193.931472 

1.054946 

23.7 

23.7 

720.338 

757.72 

1281.89 

-166.139458 

299.729320 

5.443168 

61.2 

61.2 

723.671 

317.05 

1456.00 

176.987194 

67.658507 

13.759713 

0.0 

67.7 

0.0 

1549.32 

205.968278 

176.788266 

183.822078 

26.2 

26.2 

721.833 

747.18 

1282.89 

So  version  3.3.18.0  seems  to  say  that  the  residual  veloeity  ean  be  half  the  striking 
veloeity,  roughly  one-fourth  the  striking  veloeity,  or  zero,  depending  upon  angles 
that  we  have  no  way  of  measuring  during  the  test,  even  though  the  presented  areas 
on  eaeh  of  the  orthogonal  planes  are  all  the  same. 

5.  Conclusions  

What  are  the  implieations  from  these  results?  In  many  eases,  all  we  have  are  yaw 
eard  estimates  of  the  aetual  impaet  presented  area,  so  that  at  best  we  only  have  the 
projeeted  area  on  one  plane,  the  x-y  plane.  This  may  be  adequate  for  a  eylinder,  but 
it  will  not  determine  the  orientation  of  a  euboid.  What  we  have  shown  here  is  that 
even  if  we  have  ideal  measurements  of  the  projeeted  area  on  orthogonal  planes,  this 
is  still  not  enough  to  resolve  the  orientation.  As  far  as  FATEPEN  is  eoneemed,  the 
orientation  is  speeihed  onee  the  impaet  angle  y  and  the  yaw  angle  0eff  are  known.  I 
am  not  aware  of  any  measurement  that  ean  determine  these  during  the  tests.  Thus,  I 
think  we  are  justified  in  stating  that  this  analysis  shows  that  euboids  are  not  good 
eandidates  for  EATEPEN  validation  beeause  the  results  are  highly  dependent  upon 
angles  that  we  have  no  known  way  of  measuring  and  no  way  of  uniquely  eomputing 
from  the  measured  projeeted  areas. 

What  are  the  implieations  for  trying  to  model  natural  fragments  in  EATEPEN?  If  we 
expeet  to  model  natural  fragments  as  multifaeeted  solids,  then  it  seems  we  have  a 
real  problem  sinee  we  have  no  way  of  knowing  the  many  eontaet  angles,  let  alone 
the  impaet  angle.  On  the  other  hand,  if  we  model  natural  fragments  as  the  simplest 
shape  that  has  the  same  measured  presented  area,  then  yawed  eylinders  offer  more 
promise  and  ean  be  eheeked  against  experiment. 

Part  of  the  problem  here  is  that  the  EATEPEN  eode  turns  all  shapes  into  “equivalent 

Approved  for  public  release;  distribution  is  unlimited. 


19 


cylinders”.*  So  while  we  may  be  able  to  eompute  the  eontaet  angle  of  eaeh  faee 
of  the  euboid,  that  is  not  neeessarily  the  shape  that  the  FATEPEN  eode  is  dealing 
with  at  impaet.  We  are  mueh  better  off  using  eylinders  sinee  we  ean  be  assured  that 
EATEPEN  is  using  the  same  shape.  Eurther,  it  is  more  robust  sinee  we  only  need 
to  eontrol  the  impaet  presented  area  and  there  is  only  one  eontaet  angle,  whieh  is 
the  angle  the  RCC  faee  makes  with  the  target  normal.  Onee  we  know  the  impaet 
presented  area  and  the  eylinder  dimensions,  the  impaet  angle  is  uniquely  determined. 
See  Appendix  B  for  an  implementation  of  this  proeedure. 


*See  Yatteau  et  al.,^  Section  2.1.4  on  p.  2-5  and  Section  4.1  on  p.  4-1. 
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Appendix  A.  Formula  for  an  Euler  Sequence  of  Rotations 
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Let  the  rotation  be  a  Euler  sequenee  about  3  (distinet  or  repeated)  prineipal  axes 
unit-veetors  oi,  62,  63  as  follows: 

•  First,  a  rotation  of  about  §1: 

Rl  =  -Rei(0l)-  (A-1) 

•  Seeond,  a  rotation  of  02  about  the  transformed  §2  unit  veetor: 

R2  =  Re’^{4’2)-i  (A-2) 

where  §2  =  -Rie2. 

•  Third,  a  rotation  of  03  about  the  (doubly)  transformed  03  unit  veetor: 

-R3  =  -Re"(</’3),  (A-3) 

where  §3  =  R2e'^  =  R2Ri&'i- 
The  total  eombined  rotation  is 

R  =  R3R2R1  =  -Re^'(‘/’3)-Re^(02)-Rei  (0l),  (A-4) 

where  the  rotations  are  applied  suecessively  from  right  to  left. 

Now,  the  rotation  about  §2  can  be  obtained  by  undoing  the  effect  of  the  hrst  rotation, 
performing  the  rotation  about  §2,  and  then  rotating  back: 

R2  =  Rw^{4>2) 

=  /?ii?e2(02)-Rl  ^ 

=  i?ei(0l)i?e2(02)^e/(0l)-  (A-5) 
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Similarly, 


Rs  =  -Re"  (03) 

=  -R2-Re^(03)-R2  ^ 

=  i?2-Ri-Re3(03)-Rr^-R0^  using  Eq.  A-5) 

=  [i?e,(0l)-Re3(02)-R^/(0l)]i?e,(0l)-Re3(03)i?^/(0l)[i?e,(0l)-Re3(02)-Rg/(0l)]-' 
=  i?ei(0l)i?e3(02)i?e"/(0l)^e3(0l)i?e3(03)i?r/(0l)^%(0l)^s/(02)^r/(0l) 

=  i?ei(0l)i?e3(02)i?e3(03)i?g,'(02)i?e^  (0l)-  (A-6) 

The  combined  rotation  collapses  into  something  very  simple: 


R  —  /23-R2-R1 

=  i?ei(0l)-Re2(02)-Re3(03)-Re2^(02)-R6/(0l)-Rei(0l)-Re2(02)-Re/(0l)-Rei(0l) 

=  i?ei(0l)-Re2(02)-Re3(03)-  (A-7) 


Thus,  we  get  the  result  that  the  combined  rotation  is  equal  to  the  successive  rotations 
about  the  original  unit  vectors  but  applied  in  reverse  order.  Explicitly,  this  means 


Re"(03)^e^(02)-Rei(0l)  =  RMl)  ^^2)  R&sih) 


(A-8) 


What  we  have  shown  so  far  is  more  or  less  a  plausibility  argument.  Here  is  a 
derivation  using  quaternions.  We  use  the  notation 


H'  ....  H' 

gu(0)  =  cos  -  +  usm  - 


(A-9) 


for  the  unit  quaternion  that  represents  a  counterclockwise  rotation  of  <p  radians  about 
the  unit  vector  u.  Then,  referring  to  Eqs.  A-1,  A-2,  and  A-3, 


Ri  =  gei(0i)- 


(A-10) 


where 


Y2  /  Y2 

R2  =  qe'^m)  =  cos  y  +  §2  sin  y, 


^2  =  gei(0l)e2gg,  (0l), 


(A-11) 


(A-12) 
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so  that 


and 


where 


02  /  02 
ge'(02)  =  COSy  +e2Siny 


=  COS  y  +  gei(0i)e2ggj  (0i)  sin  y 
=  gei(0i)  (^COS  y  +  §2  sin  y^  ge/(0i) 
=  gei(0l)ge2(02)ge'/(0l), 


D  /'A  ^  03  I  •  03 

^3  =  ge''(03)  =  COS  y  +  63  Sin  y , 


(A-13) 


(A-14) 


eg  =  ge' (02)e3gr0(02) 

=  ge'(02)gei(0l)e3gij^(0l)gg/(02) 

=  [gei(0l)ge2(02)g6i^(0l)]gei(0l)e3ge/(0l)[gei(0l)ge2(02)g6i^(0l)]“^ 

=  gei(0l)ge2(02)ge'/(0l)gei(0l)e3ge;^(0l)gei(0l)ge'2^(02)ge;^(0l) 

=  gei(0l)ge2(02)e3gi2^(02)g6i^(0l),  (A-15) 


so  that 


/  i  N  S^3  ,  -//  .  S^3 

ge''(03)  =  COS  y  +  63  Sin  y 

=  COS  y  +  gei(0l)ge2  (02)63^62^  (02)^6/ (0l)  sin  y 


=  gei(0l)ge2(02)  (^COS  y  +  Og  sin  y  J  {(j)2)q^,  (0l) 

=  gei(0l)ge2(02)ge3(03)ge"2^(02)ge'/(0l)- 


(A-16) 


Therefore,  the  total  eombined  rotation  is 


R  =  R3R2R1  =  geg(03)gej,(02)gei(0l) 

=  [gei(0l)ge2(02)ge3(03)ge2^(02)ge'i^(0l)][gei(0l)ge2(02)ge'/(0l)]gei(0l) 

=  gei(0l)ge2(02)ge3(03),  (A-17) 

as  was  to  be  shown. 
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The  program  in  Listing  A-1  provides  a  way  of  eheeking  this  result. 


Listing  A-1.  reverse.cpp 

1  //  reverse.cpp:  program  to  check  that  a  pitch -yaw- roll  rotation  sequence 


2 

// 

about  transformed  axes  is  equivalent  to  the  reverse  sequence 

about  fixed  axes 

3 

4 

#include  "Rotation. h" 

5 

#include  <iostream> 

6 

#include  <cstdlib> 

7 

#include  <cmath> 

8 

#include  <iomanip> 

9 

using  namespace  std; 

10 

11 

int  main(  int  argc,  char*  argv[]  )  { 

12 

13 

va::Vector  i(  1.,  0.,  0.  ),  j(  0.,  1.,  0.  ).  k(  0.,  0.,  1.  ); 

//  three  unit  vectors 

14 

va::Vector  il,  jl,  kl,  i2,  j2,  k2,  i3,  j3,  k3; 

//  transformed  unit  vectors 

15 

va::Rotation  Rp,  Ry,  Rr; 

//  rotations  for  pitch,  yaw  and  roll 

16 

double  p  =  0.,  y  =  0.,  r  =  0.; 

17 

if  (  argc  ==  4  )  { 

18 

19 

p  =  atof(  argv[l]  )  *  va::D2R; 

20 

y  =  atof(  argv[2]  )  *  va::D2R; 

21 

r  =  atof(  argv[3]  )  *  va::D2R; 

22 

} 

23 

24 

//  first,  perform  pitch  about  x-axis 

25 

Rp  =  va: : Rotation (  i,  p  ); 

26 

il  =  Rp  *  i; 

27 

j 1  =  Rp  *  j : 

28 

kl  =  Rp  *  k; 

29 

30 

//  second,  perform  yaw  about  transformed  y-axis 

31 

Ry  =  va: : Rotation (  jl,  y  ); 

32 

i2  =  Ry  *  il; 

33 

j2  =  Ry  *  jl; 

34 

k2  =  Ry  *  kl; 

35 

36 

//  third,  perform  roll  about  doubly  transformed  z-axis 

37 

Rr  =  va: : Rotation (  k2,  r  ); 

38 

i3  =  Rr  *  i2; 

39 

j3  =  Rr  *  j2; 

40 

k3  =  Rr  *  k2; 

41 

42 

cout  «  "First  done  the  conventional  way:"  «  endl; 

43 

cout  «  setprecision (6)  «  fixed; 

44 

cout  «  "i3  =  "  «  i3  «  endl; 

45 

cout  «  "j3  =  "  «  j3  «  endl; 

46 

cout  «  "k3  =  "  «  k3  «  endl  «  endl; 

47 

48 

//  perform  rotation  sequence  in  reverse  order  about  fixed  axes 

49 

va:: Rotation  R  =  va :: Rotation (  i,  p  )  *  va :: Rotation (  j,  y  ) 

*  va: : Rotation (  k,  r  ) ; 

50 

51 

cout  «  "Now,  the  same  rotations  done  in  reverse  order  about 

fixed  axes:"  «  endl; 

52 

cout  «  "i  =  '■  «  R  *  i  «  endl; 

53 

cout  «  "j  =  "  «  R  *  j  «  endl; 

54 

cout  «  "j  =  '■  «  R  *  k  «  endl  «  endl; 

55 

56 

va::Vector  v  =  3.4  *  i  -  5.7  *  j  -i-  2.3  *  k;  //  an  arbitrary 

vector 

57 

58 

cout  «  "original  vector:"  «  endl; 

59 

cout  «  V  «  endl  «  endl; 

60 

61 

cout  «  "transformed  vector  using  the  conventional  procedure: 

"  «  endl ; 

62 

cout  «  Rr  *  Ry  *  Rp  *  V  «  endl  «  endl; 

63 

64 

cout  «  "transformed  vector  using  the  reversed  order  about  fixed  axes:"  «  endl; 

65 

cout  «  R  *  V  «  endl; 

66 

67 

return  0; 

68 

} 
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For  example,  the  eommands 


C++  -02  -Wall  -std=c++ll  -o  reverse  reverse. cpp  -Im 
./reverse  65.  -137.  54. 


will  print  the  results 


First  done  the  conventional  way: 
i3  =  -0.429879  -0.021405  0.902633 
j3  =  0.591678  0.748463  0.299535 
k3  =  -0.681998  0.662832  -0.309083 

Now,  the  same  rotations  done  in  reverse  order  about  fixed  axes: 
i  =  -0.429879  -0.021405  0.902633 
j  =  0.591678  0.748463  0.299535 
j  =  -0.681998  0.662832  -0.309083 

original  vector: 

3.400000  -5.700000  2.300000 

transformed  vector  using  the  conventional  procedure: 

-6.402747  -2.814501  0.650707 

transformed  vector  using  the  reversed  order  about  fixed  axes: 
-6.402747  -2.814501  0.650707 
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Appendix  B.  Yaw  Angle  of  a  Cylinder  as  a  Function  of  Shape  Factor 
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The  dimensionless  shape  factor  7  is  defined  by  the  equation 


A,  = 


(B-1) 


where  Ap  is  the  projected  area  and  V  is  the  volume.  The  formula  for  the  dimension¬ 
less  shape  factor  of  a  right-circular  cylinder  (RCC)  as  a  function  of  yaw  angle  fy 
is 


where 


TT  L\  L 


D 


fyi 

(B-2) 

TT  L  \  TT 

(B-3) 

115)  4' 

The  yaw  angle  that  gives  the  maximum  shape  factor  is  obtained  by  setting  the 
derivative  with  respect  to  fy  equal  to  zero  and  solving  for  fy-. 


d'y 


dfy 


=  a  cos  (f)y  —  b  sin  fy  =  0, 


which  gives 


y  /  7=7max 


7=7max 


=  tan 


-1 


To  simplify  the  notation,  let  fy  denote  this  angle:  fy  =  Then, 

a  =  7max  sin  fy  and  b  =  cos  $y, 
so  that  Eq.  B-2  can  be  written  as 


Tmax  COS(0y  -  fy)  if  (j)y  >  fy 

Tmax  COS(0y  -  fy)  if  (j)y  <  fy 


where  fy  =  cos“^(6/7max)-  Solving  for  the  yaw  angle,  we  get 


I  cos  ^(6/7max) 

-fcos  ^(7/7max) 

if  7  <  6 

A=i 

I^COS  ^(6/7max) 

-  COS"^  (7/7max) 

IV 

(B-4) 


(B-5) 


(B-6) 


(B-7) 


(B-8) 


which  shows  that  we  can  easily  get  the  orientation  (yaw  angle)  of  an  RCC  from  just 
the  impact  presented  area. 
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The  code  in  Listing  B-1  implements  this  procedure  by  sampling  dimensionless  shape 
factors  from  a  lognormal  distribution  to  generate  the  appropriate  yawed  RCC.  The 
program  may  be  compiled  and  run  with  the  following  commands: 


1  C++  -02  -Wall  -std=c++ll  -o  cyl  cyl.cpp  -Im 

2  ./cyl 


Listing  B-1.  cyl.cpp 


1  //  cyl.cpp:  ImplGiTientation  of  an  algorithm  for  generating  FATEPEN  RCCs  to  represent  a  specified  shape  factor. 

2  //  Given  a  dimensionless  shape  factor,  generates  the  L/D  and  yaw  angle  for  the  RCC  to  represent  it. 

3  //  The  RCCs  are  disk-like  in  that  the  L/D  <=  Pi/4,  which  is  necessary  since  [0.5, 4. 5] 

4  //  is  the  range  of  shape  factors  for  artillery  and  rod-like  RCCs  don't  span  this  range. 

5  //  R.  Saucier,  October  2011 

6 

7  #include  <iostream> 

8  #include  <iomanip> 

9  #include  <cstdlib> 

10  #include  <cmath> 

11  #include  <random> 

12  #include  <chrono> 

13  using  namespace  std; 

14 

15  int  main(  int  argc,  char*  argv[]  )  { 

16 

17  const  int  N  =  1000;  //  number  of  samples 

18  const  double  R2D  =  180.  /  M_PI:  //  to  convert  from  radians  to  degrees 

19  const  double  SF_MIN  =  0.5;  //  minimum  shape  factor  (found  from  artillery  fragments) 

20  const  double  SF_MAX  =  4.5;  //  maximum  shape  factor  (found  from  artillery  fragments) 

21 

22  //  default  values  for  the  shape  factor  lognormal  distribution  from  122mm,  152mm  and  155mm  artillery 

23  double  mu  =  0.590494;  //  these  two  parameters  characterize  the  lognormal  shape  factor  distribution 

24  double  sigma  =  0.323433;  //  with  mode  =  1.63,  median  =  1.80  and  mean  =  1.90 

25 

26  //  ability  to  override  the  default  values  from  the  command  line  by  specifying  min  and  max  shape  factor 

27  //  that  is  meant  to  capture  95%  of  the  complete  distribution  (from  0.025  to  0.975) 

28  if  (  argc  ==  3  )  { 

29 

30  double  sfmin  =  atof(  argv[l]  ); 

31  double  sfmax  =  atof(  argv[2]  ); 

32  mu  =  0.5  *  log(  sfmin  *  sfmax  ); 

33  sigma  =  log(  sfmax  /  sfmin  )  /  (  2.  *  M_SQRT2  *  1.3859  ); 

34  } 

35 

36  unsigned  int  seed  =  std :: chrono :: high_resolution_clock :: now( ). time_since_epoch (). count () ; 

37  std::mtl9937  rng(  seed  );  //  Mersenne  Twister  engine 

38  std : : lognormal_distribution<double>  lognormal(  mu,  sigma  );  //  lognormal  shape  factor  distribution 

39  std : : uniform_real_distribution<double>  uniform(  0.069,  M_PI_4  );  //  uniform  distribution  for  L/D 

40 

41  double  a,  b,  c,  th,  sf_min,  sf_max,  sf,  l_d,  1,  d,  yaw,  V  =  1.;  //  here  we  use  a  fragment  with  unit  volume 

42  std  : :  cout  «  std  : :  setprecision  (6)  «std::  fixed; 

43 

44  for  (  int  n  =  0;  n  <  N;  n++  )  { 

45 

46  //  normally,  the  shape  factor  would  be  provided,  but  here  we  get  a  shape  factor  within  bounds  [SF_MIN,  SF_MAX] 

47  do  {  sf  =  lognormal(  rng  );  }  while  (  sf  <  SF_MIN  | |  sf  >  SF_MAX  ); 

48 

49  //  now  we  want  to  realize  this  shape  factor  with  a  cylinder 

50  do  { 

51  l_d  =  uniform(  rng  ); 

52  c  =  pow(  M_PI_4  ♦  l_d,  -2./3.  ); 

53  a  =  c  *  l_d: 

54  b  =  c  *  ^LPI_4; 

55  sf_min  =  a; 

56  sf_max  =  sqrt(  a  *  a  +  b  +  b  ); 

57  }  while  (  sf  <  sf_min  | |  sf  >  sf_max  ) ; 

58 

59  if  (  sf  <  b  )  th  =  acos(  b  /  sf_max  )  +  acos(  sf  /  sf_max  ); 

60  else  th  =  acos(  b  /  sf_max  )  -  acos(  sf  /  sf_max  ); 

61 

62  d  =  pow(  V  /  (  M_PI_4  +  l_d  ),  1./3.  ); 

63  1  =  d  +  l_d; 

64  yaw  =  th  *  R2D; 

65 

66  std:: cout  «  l_d  «  ''\t"  «  yaw  «  std::endl; 

67  } 

68 

69  return  EXIT-SUCCESS; 

70  } 
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List  of  Symbols,  Abbreviations,  and  Acronyms 

TERMS: 

FATEPEN:  Fast  Air  Target  Encounter  Penetration 
RCC:  right-circular  cylinder 
RPP:  rectangular  parallelepiped 
MATHEMATICAE  SYMBOES: 

Re{d)'.  Rotation  about  the  unit  vector  e  through  the  angle  6 
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